1. Field of the Invention
The present invention generally relates to a system and method for molecular dynamic simulation, and more particularly, to a system and method for molecular dynamic simulation which may use an algorithm having an electrostatic interaction calculating function (e.g., Ewald, particle mesh Ewald (PME), etc.) and a multiple time step function (e.g., reversible reference system propagator (RESPA, etc.)) to propagate a protein system from one set of coordinates to another in conformation space.
2. Description of the Related Art
Molecular dynamics simulations of all-atom models of proteins in water are of great current interest (J. A. McCammon and S. C. Harvey, Dynamics of Proteins and Nucleic Acids (Cambridge University Press, Cambridge, 1987); K. M. Merz, Jr. and S. M. LeGrand, The Protein Folding Problem and Tertiary Prediction (Birkhauser, Boston, 1994); J. Gunn and R. A. Friesner, Annu. Rev. Biophys. Biomol. Struct. 25. 315 (1996)). To elucidate protein folding pathways, for example, it will be necessary to simulate trajectories of duration longer than 1 microsecond (K. M. Merz, Jr. and S. M. LeGrand, The Protein Folding Problem and Tertiary Structure Prediction (Birkhauser, Boston, 1994)).
The long-range electrostatic forces in biomolecular systems make such simulations computationally intensive and lead to computational bottlenecks. Further, the dynamics in such systems usually have multiple time scales. The fast motions typical of the vibrations of intramolecular bonds require small (e.g., sub-femtosecond) time steps for stable integration of the equations of motion, and after propagating each short time step, all of the forces, including the long-range forces, must be recalculated. Since the calculation of the long-range forces is the most CPU intensive part of molecular dynamics, minimizing this part of the computational effort can lead to a great reduction of the computational cost. For this reason, considerable effort has been expended in devising methods for reducing the frequency with which the long-range forces must be recalculated, and reducing the computational effort required to compute these forces themselves.
One way to reduce the large computational cost associated with all-atom simulations is to use implicit solvent models, such as the generalized Born (GB) model, (W. C. Still, A. Tempczyk, R. C. Hawley, and T. Hendrickson, J. Am. Chem. Soc. 112, 6127 (1990)), together with stochastic dynamics with terms representing solvent friction. These models often generate very useful and interesting insights for protein folding, but it is generally difficult, if not impossible, to generate trustworthy chemical kinetic or transport information from this approach.
For realistic simulations of protein folding dynamics, all atom-based models with explicit solvent will be required. One method for reducing the computational costs of calculating the long-range forces in all-atom models is to truncate them with a spherical or minimum image truncation. This approach is very common in the literature. Unfortunately, spherical truncation or minimum image truncation is known to give rise to unphysical effects (M. Belhadj, H. E. Alper, and R. M. Levy, Chem. Phys. Lett, 179, 13 (1991); D. M. York, T. A. Darden, and L. G. Pedersen, J. Chem. Phys. 99, 8345 (1993)), and it is now widely recognized that one should not truncate the long-range electrostatic forces.
The conventional Ewald method for calculating the full Coulomb interaction in periodic systems without truncation has computational complexity that is at least of order O(N3/2). To speed up the calculation of these periodic long-range electrostatic forces, two general classes of more efficient algorithms have been developed: one is the fast multipole methods (FMM), first proposed by Greengard et al. (L. Greengard, The Rapid Evaluation of Potential Fields in Particle Systems (MIT Press, Cambridge, Mass., 1998)) and then extended to periodic systems by Francisco et al. (F. Figueirido, R. Zhou, R. Levy, and B. J. Berne, J. Chem. Phys. 106, 9835 (1997), which scales as O(N); and the other is the mesh-based Ewald methods (such as PME or the equivalent P3ME variants based on the original work of Hockney) which interpolate the point charges onto a mesh and then utilize the fast Fourier transform (FFT) to speed up the reciprocal space evaluation in the Ewald sum. These mesh-based methods scale as O(N log N). Despite the improved computational complexity of the FMM or mesh-based methods, the conventional Ewald method gives superior accuracy in determining the electrostatic forces and potentials.
Deserno and Holm (M. Deserno and C. Holm, J. Chem. Phys., 109, 7678 (1998)), have written a comprehensive and thoughtful analysis of various mesh-based Ewald methods. They compared the particle-particle particle-mesh Ewald method (or P3ME method) of Luty et al. (B. A. Luty, I. G. Tironi, and W. F. van Gunsteren, J. Chem. Phys. 103, 3014 (1995)), with the particle mesh Ewald (PME) method of Darden and with the smoothed particle mesh Ewald (or SPME) method of Essmann et al. (U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Pedersen, J. Chem. Phys. 103, 8577 (1995)), and concluded that both P3ME and SPME are considerably more accurate than PME method, and that the P3ME method is slightly more accurate than the SPME method for the same number of mesh points.
To reduce the frequency with which the long-range electrostatic forces are calculated, the reversible RESPA (r-RESPA) method (M. Tuckerman, B. J. Berne, and G. J. Martyna, J. Chem. Phys. 97, 1990 (1992)) proves very helpful. Further, in a previous publication (S. J. Stuart, R. Zhou, and Bruce J. Berne, J. Chem. Phys. 105, 1426 (1996)), a method for combining Ewald with RESPA was described. However, in that paper it is shown that a subdivision of the force, in which the real-space part of the force was included in the inner loop of r-RESPA, whereas the Fourier space part was combined into the outer loops, was inefficient. This followed because the Fourier space part contains contributions to the force which vary rapidly in time. These fast parts of the Fourier space contribution restrict the choice of the outer loop time step and thus lead to more frequent calculations of the Fourier space part than necessary.
Unfortunately, the inefficient Ewald r-RESPA strategy was adopted by others. For example, this strategy was used by Procacci et al. to combine Ewald with RESPA (P. Procacci and M. Marchi, J. Chem. Phys. 104, 3003 (1996)), and later to combine particle mesh Ewald (PME) with RESPA (P. Procacci, T. Darden, and M. Marchi, J. Phys. Chem. 100, 10464 (1996)). In this strategies, real space/reciprocal space decomposition was still used for the electrostatic forces. The Fourier, or reciprocal space (k-space), forces in PME were put in the middle of a three-level (near, medium, and long) distance-based real-space force decomposition, which leads to the real long-ranged k-space contributions being updated too often (even more often than the “long”-range real-space forces), and the short-ranged k-space contributions being updated not often enough.